﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using OSGeo.GDAL;
using OSGeo.OSR;
using System.IO;

namespace HATerrain
{
    public class TerrainFileReProjecter
    {
        /// <summary>
        /// 投影总控制函数，由于涉及参数问题未完善
        /// </summary>
        /// <param name="srcFilePath"></param>
        /// <param name="targetFileFormat"></param>
        /// <returns>重投影后文件名</returns>
        public static string Reproject(string srcFilePath,string targetFileFormat)
        {
            return null;
        }


        public static bool ReProjectToUTM(string srcFilePath,string targetFilePath,int UTMzone)
        {
            Driver driver;

            if(!File.Exists(srcFilePath)) return false;

            SpatialReference srcSR;
            SpatialReference targetSR;
            try
            { 
                Gdal.AllRegister();
                Dataset dataSet = Gdal.Open(srcFilePath, Access.GA_ReadOnly);

                TerrainHeader th;
                th = TerrainFileHelper.ReadHeader(dataSet);

                Band band = dataSet.GetRasterBand(1);

                string srcwkt = dataSet.GetProjection();
                srcSR = new SpatialReference(srcwkt);


                targetSR = new SpatialReference("");
                if (th.CoordinateNW.Y > 0)
                {
                    targetSR.SetUTM(UTMzone, 1);
                    targetSR.SetProjCS("UTM " + UTMzone + " (WGS84) in northern hemisphere.");
                }
                else
                {
                    targetSR.SetUTM(UTMzone, 0);
                    targetSR.SetProjCS("UTM " + UTMzone + " (WGS84) in southern hemisphere.");
                }
                targetSR.SetWellKnownGeogCS("WGS84");


                string targerwkt;
                targetSR.ExportToWkt(out targerwkt);


                //ss.SetProjCS("UTM 50 (WGS84) in northern hemisphere.");
                CoordinateTransformation trans = new CoordinateTransformation(srcSR, targetSR);


                double utmx, utmy;

                //经纬度转换为utm坐标
                //西南角坐标
                utmx = th.CoordinateSW.X;
                utmy = th.CoordinateSW.Y;

                double[] CoorOut = new double[3];

                trans.TransformPoint(CoorOut, utmx, utmy, 0);
                HAPoint ptSW = new HAPoint(CoorOut[0], CoorOut[1], 0);


                //西北角坐标
                utmx = th.CoordinateNW.X;
                utmy = th.CoordinateNW.Y;

                trans.TransformPoint(CoorOut, utmx, utmy, 0);
                HAPoint ptNW = new HAPoint(CoorOut[0], CoorOut[1], 0);


                //东北角坐标
                utmx = th.CoordinateNE.X;
                utmy = th.CoordinateNE.Y;

                trans.TransformPoint(CoorOut, utmx, utmy, 0);
                HAPoint ptNE = new HAPoint(CoorOut[0], CoorOut[1], 0);



                //东南角坐标
                utmx = th.CoordinateSE.X;
                utmy = th.CoordinateSE.Y;

                trans.TransformPoint(CoorOut, utmx, utmy, 0);
                HAPoint ptSE = new HAPoint(CoorOut[0], CoorOut[1], 0);


                double utmMinX = 0, utmMaxX = 0, utmMinY = 0, utmMaxY = 0;

                utmMinX = Math.Min(ptNW.X, ptSW.X);
                utmMaxX = Math.Max(ptNE.X, ptSE.X);
                utmMinY = Math.Min(ptSW.Y, ptSE.Y);
                utmMaxY = Math.Max(ptNW.Y, ptNE.Y);

                //求总长
                double x = utmMaxX - utmMinX;
                double y = utmMaxY - utmMinY;

                //求网格数
                int xSize = Convert.ToInt32(Math.Ceiling(x / 74));
                int ySize = Convert.ToInt32(Math.Ceiling(y / 75.4));

                //求间距
                x /= xSize;
                y /= ySize;

                driver = Gdal.GetDriverByName("GTiff");
                Dataset newDataset = driver.Create(targetFilePath, xSize, ySize, 1, band.DataType, new string[] { "COMPRESS=PACKBITS", "INTERLEAVE=PIXEL" });


                double[] geoInf = new double[6];
                dataSet.GetGeoTransform(geoInf);

                //求左上角坐标，转换的基准点：相对来说最合适，图像太大会有些许误差1-2/1000
                //因左上固定所以图像地面最好偏转以适应图像，否则图像信息不完整。
                trans.TransformPoint(CoorOut, geoInf[0], geoInf[3], 0);

                geoInf[0] = CoorOut[0];
                geoInf[3] = CoorOut[1];
                //设置偏转
                geoInf[2] = Math.Atan((ptNE.Y - ptNW.Y) / (ptNE.X - ptNW.X)) / Math.PI * 180;
                geoInf[4] = -Math.Atan((ptSW.X - ptNW.X) / (ptSW.Y - ptNW.Y)) / Math.PI * 180;
                geoInf[1] = x;
                geoInf[5] = -y;


                newDataset.SetProjection(targerwkt);
                newDataset.SetGeoTransform(geoInf);


                Gdal.ReprojectImage(dataSet, newDataset, srcwkt, targerwkt, ResampleAlg.GRA_NearestNeighbour, double.PositiveInfinity, 100, null, null);
                newDataset.FlushCache();
            }
            catch
            {
                return false;
            }

            return true;
        }


    }
}
